Hot-pressing process modeling for medium density 

fiberboard (MDF) 



Norberto Nigro and Mario Storti 
Centro Internacional de Metodos Computacionales en Ingenieria 
http : / /www . cimec . com . ar, <mstorti@intec . unl . edu . ar> 

February 1, 2008 
Abstract 

In this paper we present a numerical solution for the mathematical 
modeling of the hot-pressing process applied to medium density fiber- 
board. The model is based in the work of Humphrey [1982], Humphrey 
and Bolton [1989] and Carvalho and Costa [1998], with some modifica- 
tions and extensions in order to take into account mainly the convec- 
tive effects on the phase change term and also a conservative numerical 
treatment of the resulting system of partial differential equations. 

Contents 



1 Hot-pressing mathematical model 2 

1.1 Multiphase model 3 

1.2 Energy balance 3 

1.3 Steam mass balance 4 

1.4 Gas mixture mass balance 4 

2 Summary of equations and boundary conditions 5 

3 Numerical method 6 

4 Physical and transport properties 8 

4.1 Thermal conductivity 8 

4.1.1 Heat flux direction correction 8 

4.2 Permeability 8 

4.2.1 Variation of vertical permeability with board 
material density 9 

4.2.2 Horizontal permeability 9 

4.3 Steam in air diffusivity 10 

4.4 Vapor Density 10 

4.5 Saturated vapor pressure 11 

4.6 Latent heat of evaporation and heat of wetting 11 

4.7 Specific heat of mattress material 11 



4.8 Porosity 



11 



5 Numerical results 

5.1 Original numerical experiment 

5.2 Further numerical experiments 



12 

13 
13 



6 Conclusion 



18 



A List of symbols 

A.l Physical constants and quantities 

A. 2 Indices 

A. 3 Mathematical symbols 



19 

19 
20 
21 



B Derivation of the averaged energy balance equation 21 

Classification 

65M60 Finite elements, Rayleigh-Ritz and Galcrkin methods, finite 
methods 

82C70 Transport processes 

76S05 Flows in porous media; filtration; seepage 
76T30 Three or more component flows 

1 Hot-pressing mathematical model 

Hot pressing is the process in which a mattress composed of wood 
fibers and resin is cured by applying heat and pressure in a press (see 
figure 1). Continuum and batch presses do exist, and one of the main 
issues in reducing the cost of the final product is to reduce the press 
cycle time. In order to improve the heat transfer between the press 
platens and the inner layers, some amount of water is added to the 
mat. Another issue is to adjust the parameters and the temperature 
history of the cycle in order to obtain a given density profile in the 
board. Normally, it is desirable to have lower densities at the center 
of the board in order to increase the mechanical rigidity for a given 
total mass per unit area. Predicting the influence of these parameters, 
namely water content, press cycle duration and history (pressure and 
temperature) is one of the main concerns of numerical models. 

Many numerical models have been reported to help in predicting 
the influence of the process parameters in the final product. Among the 
most complete, we can find the finite difference 2D (axisymmetrical) 
model presented by Humphrey [1982] and, more recently the 3D model 
of Carvalho and Costa [1998]. Both consider conduction, phase change 
of water from the adsorbed to the vapor state and convection. The 
stress development and the determination of the density profile are 
not included in these models. 

In this paper we present a numerical model which includes all these 
features and makes some correction to the energy balance equation, as 
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presented in Carvalho and Costa [1998]. The model is based on the 
finite element model, so that it allows for a more versatile definition 
of geometry, dimensionality and, eventually, coupling with other pack- 
ages. It also will allow the use of adaptive refinement, which may be 
an important issue at the lateral borders, where the hot steam flows 
from the board to the ambient. However, this issue is not considered 
in this paper. 
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Figure 1: Hot pressing process description 



1.1 Multiphase model 

In order to avoid modeling the material down to the scale of the mi- 
crostructure (the fibers in this case), non homogeneous materials are 
solved via "averaged equations" so that the intricate microstructure re- 
sults in a continuum with averaged properties. The averaged equations 
and properties can be deduced in a rigorous way through the theory 
of mixtures and averaging operators [Whitaker, 1980]. 

1.2 Energy balance 

We will not enter in the details of all the derivations but only for the 
averaged energy balance equation, which can be found in Appendix B. 
The referred equation is 

PsC P ^ = V • (fcVT - Pv V g (C pv T + A + Qi)) - m(A + Qi), (1) 

where T is temperature, t time, C p specific heat, k thermal conduc- 
tivity, V the gradient operator, p s the density of the dry board (solid 
phase), p v vapor density, V g the volume averaged gas velocity, i.e. 

V 3 = ev g (2) 

where v g is the velocity averaged on the phase, see Appendix B), C pv 
specific heat of vapor, m evaporation rate, A latent heat of vaporization 
of free water, and Qi adsorption heat. The main difference between 
this equation and that presented by both Carvalho and Costa [1998] 
and Humphrey [1982] is in the addition of the water evaporation heat 
term in the convection term instead of considering the phase change 
effect only on the temporal term. This term should be included because 
both phases, the solid material and the vapor are in relative motion 
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and we think that its influence is not negligible in a high temperature 
process. 



1.3 Steam mass balance 

Carvalho and Costa [1998] proposed the following steam mass conser- 
vation equation 



. MM W 
m= _ cV . 



• |V,^ 



(3) 



where MM W is the molecular weight of water, R the gas constant, e 
board porosity, D v diffusivity of water vapor in the air/vapor mixture 
and P v vapor partial pressure. Considering that the steam is treated 
as an ideal gas, then 

Y = MM W _1 %, (4) 

so it may be written, assuming e constant, as 

m = V • [-eD v Vp v + V g p v ] . (5) 

This expression is preferable to (3) because it is written in a conserva- 
tive form that is more agreeable for a numerical treatment. The left 
hand side term represents the mass interfacial transport and those in 
the right hand side take into account the mass diffusion and the mass 
convection. However, it should be noted that this last expression does 
not have a temporal term as every consistent balance equation does. 
For example, if evaporation is not considered, then (5) is valid only for 
a steady situation, which is not in general the case. Then, we rewrite 
the steam mass balance as 

4^ = V [eD v Vp v - V gPv ] + rh. (6) 

This is another difference between our model and that proposed by 
Carvalho and Costa [1998]. 

1.4 Gas mixture mass balance 

Finally, because the gas phase is composed of two main constituents, 
steam and air, we may use an additional equation for the mass trans- 
port of the whole gas phase. Carvalho and Costa [1998] considered 

dP I f K g P \ rh P dT 

where P is the pressure of the gas phase and K g the board permeability 
tensor. Again, assuming ideal gas law as the state equation for this 
phase, 

d (P\ R ( K q \ mR 
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we arrive to 

e^= -V-^V.J + m. (9) 



In order to close the system of equations we need to introduce 
a relationship between rh and (dP v /dt). Consider the steam mass 
balance (6) and the relation 

p s U = ep v + p L , (10) 

that represents the fact that the board moisture content H is composed 
of vapor and bound water p^. If we assume that no liquid phase is 
considered, then bound water may be transferred to the gas phase only 
(solid to steam). So 

(id 

and then 

Ps at " 6 dt dt (12) 

= V • [eD v V Pv - V gPv ] . 

The air mass balance equation can be obtained by substracting (6) 
from (9) 



= V • [-eD v V Pv - V gPa ] . (13) 



Due to the fact that the mean macroscopic diffusive fluxes should be 
null 

D v V Pv + D a V Pa = 0, (14) 

the air mass balance equation is transformed in the following expression 

dp a 
dt 



^ = V ' [eDaVpa - V gPa ] , (15) 



which is very similar to (6) but here valid for the air. Obviously, the 
air transport equation has no evaporation term. 

2 Summary of equations and boundary 
conditions 

In order to clarify the mathematical model that is finally used for 
the simulation of hot-pressing process we present the following brief 
summary of partial differential equations. 



Energy balance equation 
dT 

PsC P -^ = V • (fcVT - Pv V g (C pv T + A + Qi)) - m(X + Q t ). 

(16) 
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• Water content balance equation 



<9H 

p s -^ = V-[eD v Vp v -V g p v ]. (17) 



• Air mass balance equation 

6 ^t = V ' [eDaVf>a ~ Y 3Pa] ■ (18) 

The boundary conditions are the following: 

• At the press platen: 

T = Tpi at0 n(i), air/ water mixture in equilibrium with platen temperature 

V g ■ n = 0, no mass flow across the platen 

dp 

— — = 0, no air diffusion across the platen 
Oz 

dp 

— - = 0, no vapor diffusion across the platen. 
Oz 

(19) 

• At the center line (r = 0), axial symmetry for all variables 

dT dH dp a , , 

— =0, -^-=0, -^=0. 20 
or or Or 

• At the mid plane (z = 0), symmetry for all variables 

dT dH dp a , , 

— =0, -^-=0, -^=0. 21 

Oz oz Oz 

• At the exit boundary (r = R ex t), 

dT 

— — = 0, null diffusive heat flux 
Or 

Pv = Pv,atm, cquil. with external air/water mixture, ^ ' 
P a = Pa.atm, equil. with external air/water mixture. 

3 Numerical method 

The above system of equations contains three main unknowns, the 
temperature, the moisture content and the air density representing 
the dependent variables of the problem also called the state variable. 
In this work we have used as independent variables the time and two 
spatial coordinates (3D problems may be computed much in the same 
way) . Due to the physical and geometrical inherent complexity of this 
problem this may be computed only by numerical methods. For the 
spatial discretization we have employed finite elements with multilinear 
elements for all the unknowns. Due to the high convective effects the 
numerical scheme was stabilized with the SUPG (for "Streamline Up- 
wind - Petrov Galerkin") formulation (see Brooks and Hughes [1982]), 
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otherwise spurious oscillations arise. Once the spatial discretization is 
performed the partial differential system of equations is transformed 
into an ordinary differential system of equations like 

U = R{U), (23) 

where U is the state vector containing the three unknowns in each 
node of the whole mesh. So, the system dimension is 37V where N 
is the number of nodes in the mesh. The numerical procedure is as 
follow: Knowing the state vector at the current time (t n ), i.e. Uj(t n ) = 
[Tj, H f, p a ,j](t n ) where j represent an specific node in the mesh. To get 
the residual, right hand side of (23) the following steps should be done: 

• Obtain air pressure from the gas state equation ((T,p a ) — > P a ). 

• Obtain relative humidity from sorption isotherms 
((T,H)->HR). 

• Compute saturated vapor pressure from Kirchoff expression 
(38) (T— > P sat ). 

• Compute vapor pressure P v = HRP sat . 

• Obtain vapor density in air from vapor state equation 
(T, P v ) -» Pv 

• Compute coefficients from additional constitutive laws (P, H, T) 

* {D , Kg, k x ,y i Cp) . 

• Compute gradients of T, P from nodal values at the Gauss 
points using the finite element interpolation. 

• Assemble the element residual contributions in a global vector. 

Once the whole residual vector at t = t n is computed the unknowns 
variables at the next time step is updated with 

jjn+l =U n + Ati?([/™). (24) 

This kind of scheme, called explicit integration in time is very simple 
to be implemented but it has two major drawbacks, one is the limi- 
tation of the time step to ensure numerical stability and the other is 
the bad convergence rate for ill-conditioned system of equations. In 
this application the last disadvantage is very restrictive because the 
characteristic times of each equation are very different. To circumvent 
this drawback we have implemented an implicit numerical scheme 

U n+1 - AtR(U n+1 ) = U n , (25) 

where the residue is computed at the new state variable U(t n+1 ) in- 
stead of using the current value U(t n ). The non-linearities and the time 
dependency of the state vector makes this implementation more diffi- 
cult and time consuming but more stable. This non-linear equation in 
jjn+i j g so i vec j ^ the Newton method, which requires the computation 
of the Jacobian 



In order to avoid the explicit computation of the Jacobian, we compute 
an approximate one by finite differences. 



J S3 J, 



num — 



R{U + 5U)-R(U) 
SU 



(27) 



This can be done element-wise, so the cost involved is proportional to 
the number of elements in the mesh. In our application and despite 
of the ill-condition of the problem we have found a good convergence 
at each time step, in average 4 iterations per time step to reduce 10 
orders of magnitude in the global residue. 

4 Physical and transport properties 
4.1 Thermal conductivity 

Following Humphrey [1982] the influence of board density, moisture 
content and temperature on the thermal conductivity is considered as 
an independent correction factor obtained experimentally 



where k z is the thermal conductivity in the pressing direction in 
W/m K and p s is the oven dry density of the material in Kg/m 3 . 
The moisture correction factor F Kt n (Kollmann and Malmquist [1956] 
and Humphrey [1982]) assumes H, the moisture content of the board 
material, in %, and the temperature correction factor F K j- (Kuhlmann 
[1962] and Humphrey [1982]) assumes T in °C. 

4.1.1 Heat flux direction correction 

By far the greater part of conductive heat translation takes place in 
the vertical plane. However the energy lost from the mattress is largely 
the result of radial vapor migration from the center toward the atmo- 
sphere. The associated horizontal relative humidity gradient lead to a 
horizontal temperature gradient. Even though this gradient is always 
lower than the vertical one its influence should be taken into account 
if multidimensional analysis is required. Ward and Skaar [1963] made 
experimental measurements and they observed that at a first glance a 
factor of approximately 1.5 may be a good initial guess before doing 
some extra measurements. Then 



4.2 Permeability 

The evaporation and condensation of water changes the vapor density 
and consequently its partial pressure in the voids within the composite. 




F K , H F KyT (1.172 x 1(T 2 + 1.319 x 1(T 4 p a ) 
l + 9.77x 10" 3 (i?-12), 
(T - 20) x 1.077 x 10 3 + 1, 



(28) 




(29) 



8 



A vertical pressure gradient leads to the flow of water vapor from the 
press platens toward the central plane of the board. At the same 
time an horizontal vapor flow is set up in response to the pressure 
gradient established in the same direction. The relation between the 
pressure gradient and the flow features may be assigned to the material 
permeability. Permeability is a measure of the ease with which a fluid 
may flow through a porous medium under the influence of a given 
pressure gradient. Different mechanisms may be involved in this flow, 
a viscous laminar flow, a turbulent flow and a slip or Knudsen flow. 
In this study only the first type is included with the assumption that 
Darcy's law is obeyed. This may be written as 

V 9 = -^Vp (30) 
Ms 

where Vp is the driving force and v 9 is the flow variable. K g /fi g is 
called the superficial permeability with n g the dynamic viscosity of the 
gas phase and K g the specific permeability. 

The kinetic theory of gases suggests that at normal pressure the 
viscosity is independent of pressure and it varies as the square root of 
the absolute temperature, 

H oc Vf 10 3 < p < 10 6 . (31) 

Corrections with the absolute temperature are often considered by 
Sutherland law 

B T 3 / 2 

(32) 

with B and C characteristic constants of the gas or vapor with /i in 
Kgm~ x s _1 . Values for B and C are available from Keenan and Keyes 
[1966]. For this application the following expression was adopted, 

„ = 1.112 xl0-*x(™^, (33) 
P (T + 3211.0) v ; 

with T in °C. 



4.2.1 Variation of vertical permeability with board ma- 
terial density 

Even though we consider that the press is closed and the density profile 
is set up and fixed it is included here some conclusions from Humphrey 
[1982] with results obtained by Denisov et al. [1975] for 19 mm boards 
and others from Sokunbi [1978]. The data to be fitted are the following 



4.2.2 Horizontal permeability 

Sokunbi measures included in figure 2.7 of Humphrey [1982] shows the 
relation between the board thickness in mm with horizontal perme- 
ability. For approximately 15mm board thickness and p s — 586Kg/m 3 
the horizontal permeability is 59 times the vertical value, in agreement 
with the values assumed by Carvalho and Costa [1998]. 



9 



Mean density 
[Kg/m 3 ] 


Mean vertical 
permeability 
[m 2 x 10 15 ] 


425 


64 


475 


40 


525 


24 


575 


16 


625 


11 


675 


7 


725 


5 


775 


3 


825 


2 


875 


2 



Table 1: Table I: Vertical permeability density correction data 



4.3 Steam in air diffusivity 

The interdiffusion coefficient of steam in air can be calculated from the 
following semi-empirical equation [Stanish et al., 1986], 

^, 20xl0 - 5 (i2^)(_!_) (34) 

where the diffusivity is in m 2 /sec, pressure in N/m 2 and T in °K. 

4.4 Vapor Density 

For the pressure range likely to occur during hot pressing (between 
10 3 and 3 x 10 5 N/m 2 ) a linear relationship between saturated vapor 
pressure and vapor density may be assumed. Fitting experimental data 
Humphrey [1982] proposed the following expression 

p v = P sat 6.0 x HT 8 HR, (35) 

with p v in Kg/m 3 , P sat in N/m 2 and the relative humidity HR in %. 
This can be deduced from the relative humidity definition 

HR = P v /P S au (36) 
and applying ideal gas law for gaseous phase 

— = R/MM W T. (37) 

Pv 

Taking R = 8314J/Kmol/K and MM W = 18Kg/Kmol with T w 360°K 
we obtain (35). 
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4.5 Saturated vapor pressure 

Following the Kirchoff expression with data presented by Kecnan and 
Keyes [1966] we include here the following equation, 

logic Psat = 10.745 - (2141.0/(T + 273.15)), (38) 

with P sat in N/m 2 and T in °C. 



4.6 Latent heat of evaporation and heat of wetting 

Using Clausius-Clapeyron equation in differential form and after some 
simplification the latent heat of vaporization of free water may be writ- 
ten as 

A = 2.511 x 10 6 - 2.48 x 10 3 T, (39) 

with A in J/Kg and T in °C. 

For the differential heat of sorption we follow Humphrey [1982] that 
used (see Bramhall [1979]) 

Qi = 1.176 x 10 6 c~ 015H , (40) 

with Qi in J/Kg and H in %. 



4.7 Specific heat of mattress material 

It is computed by adding the specific heat of dry wood and that of 
water according to the material porosity and assuming full saturation. 
The specific heat of dry mattress material is taken as 1357J/Kg/K and 
the specific heat of water has been taken to be 4190J/Kg/K. From Siau 
[1984] the expression for specific heat of moist wood is 

C p = 4180°- 268 + - 0011(T - 273 - 15) + g , (41) 
1 + H 

T in °K and H in %. 



4.8 Porosity 

According to Humphrey [1982] the volume of voids within the region 
may be computed with 

e=%^ = (l-f), (42) 
V p s 

where e is the porosity, p is the density of the region and p s is the dry 
density of the board material. 

In Carvalho and Costa [1998] they included the expression from 
Suzuki and Kato [1989] 

e = l- P& l,P \l Vr,P \ (43) 

J- T Vr 
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Figure 2: Finite element mesh 



where p r is the cure resin density, pf is the oven dry fiber density and 
y r is the resin content (resin weight/board weight). In Carvalho and 
Costa [1998], the authors used y r = 8.5% with p/ = 900Kg/m 3 and 
p r = 1100Kg/m 3 . 

5 Numerical results 

In this section we present some results that can be compared with 
those reported in Humphrey [1982]. This numerical experiment allows 
the validation of the mathematical model and its numerical implemen- 
tation for future applications to hot-pressing process simulation and 
control. This experiment consists of a round fiberboard of 15 mm of 
thickness and 0.2828 m of radius that according to its axisymmetri- 
cal geometry needs as spatial coordinates only the radius r and the 
axial coordinate z assuming no variation in the circumferential coordi- 
nate. The axisymmetrical domain is discretized in 20 x 20 elements in 
each direction with a grading toward the press platen and the external 
radius as may be visualized in figure 2. In order to follow the same 
assumptions as in that work we fixed the air density to a very low value 
(p a ~ 10~ 6 ), as if the press was close with no air inside the fiberboard. 
The boundary conditions are as in section 2. 

For the press platen temperature we have applied a ramp from 30 C 
at t = to 160°C at t = 72 seconds with a least square fitting from data 
in Humphrey [1982]. As initial conditions we have assumed a uniform 
temperature of T{t = 0) = 30 C in the whole solid material with a 
uniform moisture content of H = 11%. The external atmosphere was 
considered to be at T atm = 30 C, HR a tm = 65%, in such a way that 
the internal moisture content is, at the initial state, in equilibrium with 
the external atmosphere. 

In the next section we include the results obtained by using the 
model above cited to the original problem presented in Humphrey 
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[1982]. Next we present some other experiments showing some phe- 
nomena that deserve more attention for futue studies. 

5.1 Original numerical experiment 

Figure 3 shows the temperature distribution in time and with the 
axial coordinate at r = (ccntcrline). We can note the penetra- 
tion in the axial direction of the temperature profile in time for 
t = 1,10,50,100,200,300 and 400 seconds. Figure 4 shows the same 
kind of plot for moisture content. For r not to close to the external ra- 
dius, the problem is almost one dimensional in the z direction. As the 
thermal front penetrates into the board water evaporates. This vapor 
advances to lower pressure regions near the symmetry plane and, as it 
encounters lower temperatures, it condenses releasing heat. This pro- 
cess can be clearly seen from the wave in moisture content (figure 4) 
exceeding the initial water content of 11% and results in an improve- 
ment in the heat transfer with respect to the pure conduction case. 
Also this phenomenon is responsible of the change in curvature of thee 
temperature curves, mainly at t — 10, 50, and 100 sec (see figure 3). 
The total water content in the board at a given instant can be found 
by integrating the bound water content and the water in vapor phase. 
However, this last is negligible. We can see in figure 4 that the depres- 
sion in water content near the board (for instance at t = 400sec), is 
larger that the water enrichment in near the center plane. This is due 
to water migration from the center of the board to the external radius, 
where it flows to the external atmosphere. The following figures show 
similar plots but at different locations, 

• Figure 5: temperature at external radius 

• Figure 6: moisture content at external radius 

• Figure 7: temperature at axial centerline (z = 0) 

• Figure 8: moisture content at axial centerline (z = 0) 

• Figure 9: moisture content at press platen 

Figure 10 shows several isotherms at t = 200 seconds distributed 
in the r, z plane. 

Figure 11 shows a three-dimensional view for moisture content rep- 
resented by the third coordinate axis at t = 200 seconds as a function 
of r, z. 

These results are in good agreement with those presented by 
Humphrey [1982]. However the cited author did not present his re- 
sults at some locations that in our opinion should be treated with 
some care, for example at the external radius. 

5.2 Further numerical experiments 

In Humphrey [1982], results for moisture and temperature in the ver- 
tical and radial directions at both central planes, r — and z = 
respectively are included. No mention about the vertical distribution 
at r — R = 0.2828m or about the moisture content at press platen. 
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Figure 3: Axial temperature profile at r = 
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Figure 4: Axial moisture content profile at r = 
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Figure 5: Axial temperature profile at r = R = 0.2828m 
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Figure 6: Axial moisture content profile at r = R = 0.2828m 
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Figure 8: Radial moisture content profile at z = 
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Figure 9: Radial moisture content profile at press platen 
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Figure 10: Isotherms at t=200 seconds (z and r axis not to scale.) 
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Figure 11: 3D view for moisture content at t=200 seconds 

Moreover, he had used a uniform mesh of 10 x 10 elements without 
showing what happen at the last annuli of elements corresponding to 
the external radius. 

Our results present some overshooting in the temperature profile 
very close to the radial exit contour and at the first moment we though 
about a spurious numerical problem, but it is due to large variations 
of the magnitude of vapor pressure and density at the boundary. In 
typical runs, vapor pressure varies from near 2 atm in the center of the 
board to 0.01 atm at the external radius. We think that this problem 
will be fixed if we solve for the air density also, but then a very fine 
grid will be necessary at the exit boundary, since large variations of 
the vapor molar fraction are expected, (see figure 12). Molar fraction 
varies from nearly 1 at the interior of the mat to a 2% at the external 
atmosphere. This variation is produced in a thin layer of width 5 
proportional to the diffusivity of vapor in air which is very small. 

6 Conclusion 

We presented a numerical model for the heat and mass transfer in the 
hot-pressing model of a MDF fiberboard. The model includes convec- 
tive effects on the phase change term and also a conservative numerical 
treatment of the resulting system of partial differential equations. Con- 
vective effects are responsible of an increase in heat transfer from the 
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Total pressure 




air 
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Figure 12: Boundary layer in vapor partial pressure at external radius 

platen to the center of the board due to water vapor evaporation and 
condensation. Two-dimensional simulations allow to estimate border 
effects. 
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A List of symbols 

A.l Physical constants and quantities 

D: diffusivity [m 2 /s] 
At: time step [s] 

e: material porosity (dimcnsionless) 
Fk'. correction factor for thermal conductivity (dimensionless) 
Fo: Fourier number (dimensionless) 

h: enthalpy 

H: water content in % weight 
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HR: 



J: 



k: 



relative humidity [%] 

Jacobian matrix 

thermal conductivity [W/m°K] 



K g : specific permeability [m 2 ] 

A: latent heat of vaporization of free water [J/Kg] 

m: evaporation rate [Kg/m 3 s] 

MM: air molecular weight [Kg/Kmol] 

\i: dynamic viscosity [Kg/ms] 

N: number of nodes in the finite clement mesh 

P: pressure [N/m 2 ] 

Qi: adsorption heat [J/Kg] 

R: gas constant for air (8314KJ/Kmol K 

R: vector of residuals for the discrete model 

r: radial coordinate [mm] 

p: density [Kg/m 3 ] 

T: temperature [°K] 

U: state vector for the discrete model 

V: volume [m 3 ] 

v: gas mixture velocity [m/s] 

Koids: volume of voids [m 3 ] 

y r : resin content [weight %] 

z: coordinate normal to the plate [mm] 

A. 2 Indices 

eff: effective quantities (averaged for the gas/solid mixture) 

/: fiber 

g: gas phase, (air/water mixture) 

j: nodal index 

L, I: bound water 

platen: quantity evaluated at the press platen 

r: resin 

ref: reference state 

s: solid phase 

sat: saturated atmosphere 

v: water vapor 

w: water 
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A. 3 Mathematical symbols 

V: gradient (nabla) operator 
( ): temporal derivative 

B Derivation of the averaged energy bal- 
ance equation 

The microscopic energy balance equation in the gas phase is 

^ (ph) + V • (pwh) = - V • (fcVT) . (44) 

where h is enthalpy. For the other phases (solid and bound water) a 
similar expression holds, but neglecting the advective term. Apply- 
ing the volume average operator [Whitaker, 1980], we arrive to the 
following equation averaged on the gas phase 



d_ 

dt 



-(e g (ph) 9 ) + W-(e g (phv) 9 ) = V ■ (e fl (kWT) 9 ) + Q g , (45) 



e g is the volumetric fraction of phase g (i.e. gas) and (X) 9 is the 
average of quantity X on the volume occupied by phase g 

(X) 9 = ±- [ Xdn. (46) 



Jn g 

The term Q g is the total enthalpy flux through the solid-gas interface 

r 

Q g = J{ph) g {v-w)-ndT, (47) 

where (X) g is the value of property X on the g side of the interface 
and w is the velocity of the interface. Assuming that h g is constant 
on all f2 g (for a certain volume control) then 



Q a = W 9 ^(p) 9 (v-w).ndl\ 



Ir"" ' (48) 

= (h) 9 rh, 

where rh is the rate of mass of water being evaporated. A common 
drawback of averaged equations is that, when products of variables 
like ph appear in the microscopic equation, the average of the product 
(ph) 9 is obtained in the averaged equation. Now, it is not true that 

(ph) 9 = (p) 9 (h) 9 , (49) 

so that the averaged equation contains more unknowns than the orig- 
inal equation. A common assumption is that no correlation exists 
between variables and so (49) is approximmately valid. This can be 
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justified, for instance, if the variations of each quantity around the 
mean is small. 

Then, applying the volume average operator over the gas, solid, and 
bound water phases and assuming no correlations between variables we 
obtain the following averaged equations for the three phases 

d 

Ol^gPghg) + v ' ( e gPg v g h g) = V ■ (e fl (fcVT) 9 ) - mh g gas phase, 
d 

— (e s p s h s ) — V • (e s (fcVT) s ) solid phase and , 

d 

-^(eipihi) =rhhi liquid phase. 

(50) 

Also, the average operator is dropped from here on, and a subindex g 
or s implies averaging on that phase. Also, V g the volume averaged 
gas velocity is 

1 



V g = - / vdJ] = t9 v s . (51) 



Note that in the body of the text V g is used instead of vg. Now, 
summation of these three equations gives 

d 

-g^( € gPg h g + e sPsh s + eipihi) + V • (e g v g p g h g ) = 

V(e s (fcVT) 9 + e s (fcVT) s ) + m(h t - h g ). (52) 

Now, 

e g (kVT) 9 + e s (kVT) s = k eS V (T) , (53) 

where fc e ff is the average conductivity of the solid+water+gas mixture. 
The gas is assumed as an ideal mixture, so that the enthalpy is the 
sum of the enthalpy of its constituents, and neglect the contribution 
of the air constituent so that 

p g h g = p a h a + p v h v w p v h v . (54) 

Taking a reference state for the entalphy at a point on the adsorbed 
state 

h v = C pv (T-T ref ) + X + Qi. (55) 

We also neglect the entalphy of the gas phase with respect to the 
solid+water phases and put 

e s p s h s + eipihi = {pC p ) c s (T - T re f), (56) 

where p e ff, and C Pcff are averaged properties for the moist board, as 
a function of temperature and moisture content. Finally, the averaged 
equation is 

^[(pC p )e ff (T - T ref )] + V • [e gPv v g (C pv (T - T ref ) + A + Q,)] = 

V ■(k cB T)-m(X + Q l ) (57) 
This is equivalent to (1) through relation (2) and assuming T ro f = 0°C. 
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